COMMENT ⊗   VALID 00002 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	ENTRY GETBDR,STOBDR,SHPROPS,shape,RECONSTRUCT
C00012 ENDMK
C⊗;
ENTRY GETBDR,STOBDR,SHPROPS,shape,RECONSTRUCT;
BEGIN "shape"
REQUIRE "36A" COMPILER!SWITCHES;
REQUIRE "VISLIB" SOURCE!FILE;
source!v(extitm);
external integer itemvar tbdrlen,tshapes,ilocv,jlocv;
EXTERNAL INTEGER MODIFIED;
require "⊂⊃<>" delimiters;
define notiming=⊂TRUE⊃;
source!l(saitim);
own integer inittsh;
SIMPLE PROCEDURE SHAPEINIT;
    BEGIN
    IF INITTSH THEN RETURN;
    tshapes←dcrepro("TSHAPES",9,0);
    ilocv←dcrepro("ILOCV",2,0);
    jlocv←dcrepro("JLOCV",2,0);
    tbdrlen←dcrepro("TBDRLEN",2,0);
    inittsh←-1;
    end;

simple internal integer procedure getbdr(string itemvar reg; integer bigmsk);
    begin "getbdr"
    integer itemvar x1;
    if tbdrlen ⊗ reg ≡ BIND x1
	THEN BEGIN
	    if props(x1) then modified←props(x1);
	    if bigmsk
		then return(lhalf(<datum(x1)>))
		else return(rhalf(<datum(x1)>))
	    end
	else return(0);
    end "getbdr";

simple internal procedure stobdr(string itemvar reg; integer bigmsk, bdrl);
    begin "stobdr"
    integer itemvar x1;
    if tbdrlen ⊗ reg ≡ BIND x1
	then if bigmsk
	    then datum(x1)←(bdrl LSH 18) LOR rhalf(<datum(x1)>)
	    else datum(x1)←(lhalf(<datum(x1)>) LSH 18) LOR bdrl
	else dwrite(tbdrlen,reg,dnew(if bigmsk then (bdrl LSH 18) else BDRL));
    end "stobdr";
simple internal integer procedure shprops(string itemvar reg;safe real array an,thn,bn,psin);
    begin "shprops"
    safe real array itemvar raiv;
    integer size1,i;
    shapeinit;
    if ¬(tshapes ⊗ reg ≡ BIND raiv) then return(-1);
    size1←arrinfo(datum(raiv),2);
    for i←1 thru size1 do
	begin
	an[i-1]←datum(raiv)[i,1];
	thn[i-1]←datum(raiv)[i,2];
	bn[i-1]←datum(raiv)[i,3];
	psin[i-1]←datum(raiv)[i,4];
	end;
    return(0);
    end "shprops";

INTERNAL integer PROCEDURE shape(string itemvar reg; INTEGER numharm,bigmsk; safe real array an,thn,gamn,bn,psin,deln; integer force);
	BEGINtim(shapeFOLLOW,
	<SAFE INTEGER ARRAY NEIGHBORS[0:7];
	safe real array axn,bxn,ayn,byn[0:numharm];
	INTEGER REGNUM,M,N,INITI,INITJ,P,OLDVAL,RWS,COLS,I,J,ST,TEMP,numpt,val,offi,offj,ptnum,ind,ibuf,flag;
	string itemvar mname;
	integer itemvar x1;
	string maskn;
	real pi2m,pi2mk,ttrig,tang,xs,ys;>)
	if force=0 then if shprops(reg,an,thn,bn,psin)=0 then blreturnit(-1);
	shapeinit;
	if dmask ⊗ reg ≡ BIND mname
	    then begin
		indmp(getdev(maskn←datum(mname),if bigmsk then "BIG" else "MSK"),maskn,ibuf←fndbuf,flag←-2);
		if flag then blreturnit(0);
		end
	    else blreturnit(0);
	RWS←ROWS(IBUF);  COLS←COLMS(IBUF);
	uptoval(i←1,j←1,1,ibuf);
	initi←i;  initj←j;
	ST←0;
	offi←isubst(ibuf);
	offj←jsubst(ibuf);
	if ¬(ilocv ⊗ reg ≡ BIND x1) then dwrite(ilocv,reg,dnew((offi LSH 18) LOR (rws+offi-1)));
	if ¬(jlocv ⊗ reg ≡ BIND x1) then dwrite(jlocv,reg,dnew((offj LSH 18) LOR (cols+offj-1)));
	if (numpt←getbdr(reg,bigmsk))=0
	    then begintim(first follow)
		numpt←0;
		BDRPRE(N,I,J,IBUF,RWS,COLS,NEIGHBORS);
		WHILE TRUE DO
			BEGIN "loop"
			TEMP←(N+ST) MOD 8;
			I←I+ICASEV(TEMP);
			J←J+JCASEV(TEMP);
			numpt←numpt+1;
			IF I=INITI AND J=INITJ THEN done "loop";
			BDRPOST(N,ST,TEMP,I,J,IBUF,RWS,COLS,NEIGHBORS);
			END "loop";
		stobdr(reg,bigmsk,numpt);
		endtim(first follow)
	    else ;
	! do it again to get the Fourier Xform stuff;
	begintim(Calculate terms)
	begintim(Actual terms)
	arrclr(axn);  arrclr(bxn);
	arrclr(ayn);  arrclr(byn);
	pi2m←3.1415927*2.0/numpt;
	endtim(Actual terms);
	ptnum←ST←0;
	i←initi;  j←initj;
	BDRPRE(N,I,J,IBUF,RWS,COLS,NEIGHBORS);
	WHILE TRUE DO
		BEGIN "loop2"
		TEMP←(N+ST) MOD 8;
		I←I+ICASEV(TEMP);
		J←J+JCASEV(TEMP);
! here goes the computations;
		begintim(Actual terms)
		ptnum←ptnum+1;
		pi2mk←pi2m*ptnum;
		for ind←0 thru numharm do
		    begin
		    axn[ind]←axn[ind]+i*(ttrig←cos(tang←(pi2mk*ind)));
		    ayn[ind]←ayn[ind]+j*ttrig;
		    bxn[ind]←bxn[ind]+i*(ttrig←sin(tang));
		    byn[ind]←byn[ind]+j*ttrig;
		    end;
		endtim(Actual terms);
		IF I=INITI AND J=INITJ THEN done "loop2";
		BDRPOST(N,ST,TEMP,I,J,IBUF,RWS,COLS,NEIGHBORS);
		END "loop2";
	! now get the final computations;
	begintim(Actual terms)
	gamn[0]←an[0]←axn[0]/numpt;
	deln[0]←bn[0]←ayn[0]/numpt;
	thn[0]←atan(bxn[0]/axn[0]);
	psin[0]←atan(byn[0]/ayn[0]);
	for ind←1 thru numharm do
	    begin
	    pi2m←3.1415927*ind;
	    ttrig←sin(tang←(pi2m/numpt));
	    pi2m←2.0/pi2m;
	    an[ind]←pi2m*ttrig*sqrt(axn[ind]*axn[ind]+bxn[ind]*bxn[ind]);
	    bn[ind]←pi2m*ttrig*sqrt(ayn[ind]*ayn[ind]+byn[ind]*byn[ind]);
	    thn[ind]←atan2(bxn[ind],axn[ind])-tang;
	    psin[ind]←atan2(byn[ind],ayn[ind])-tang;
	    end;
	endtim(Actual terms);
	endtim(Calculate terms);
	begin "add relations"
	safe real array rarr[1:numharm+1,1:4];
	for i←0 thru numharm do
	    begin
	    rarr[i+1,1]←an[i];
	    rarr[i+1,2]←thn[i];
	    rarr[i+1,3]←bn[i];
	    rarr[i+1,4]←psin[i];
	    end;
	dwrite(tshapes,reg,dnew(rarr));
	end "add relations";
	frebuf(ibuf);
	blreturnit(numpt);
	endtim(shapeFOLLOW);

SIMPLE internal integer procedure reconstruct(string itemvar reg; integer numharm; safe real array an,thn,gamn,bn,psin,deln);
	begintim(reconstruct,
	<integer nbuf,i,n,rws,col,numpt;
	integer itemvar x1,x2;
	real xs,ys,pi2m,tang;>)
	SHAPEINIT;
	if ¬(ilocv ⊗ reg ≡ BIND X1) then blreturnit(-1);
	if ¬(jlocv ⊗ reg ≡ BIND X2) then blreturnit(-1);
	if shprops(reg,an,thn,bn,psin) then blreturnit(-1);
	getbuf(rws←(rhalf(<datum(x1)>)-lhalf(<datum(x1)>)+1),col←(rhalf(<datum(x2)>)-lhalf(<datum(x2)>)+1),1,nbuf←fndbuf);
	putsub(lhalf(<datum(x1)>),lhalf(<datum(x2)>),nbuf);
	if (tbdrlen ⊗ reg ≡ BIND x1) then numpt←lhalf(<datum(x1)>);
	for i←1 thru numpt do
	    begin
	    xs←an[0];
	    ys←bn[0];
	    pi2m←3.1415927*2*i/numpt;
	    for n←1 thru numharm do
		begin
		xs←xs+an[n]*cos((tang←(pi2m*n))-thn[n]);
		ys←ys+bn[n]*cos(tang-psin[n]);
		end;
	    putpnt((xs MAX 1) MIN rws,(ys MAX 1) MIN col,-1,nbuf);
	    end;
	blreturnit(nbuf);
	endtim(reconstruct);
END "shape"